#!/bin/bash

ref12=/lustre/home/acct-medfzx/medfzx-lkw/project/CAH/data/refseqMy/hg19/hg19

Time=TRF-20241105-L-01-2024-11-091725

outputPath=/lustre/home/acct-medfzx/medfzx-lkw/project/m4/RESULT_BASH_final_20241102/${Time}
samplePath=/lustre/home/acct-medfzx/medfzx-lkw/project/m4/data/fastq/${Time}
BED=/lustre/home/acct-medfzx/medfzx-lkw/project/m4/bed/m4_2.bed  # 目标区域文件

mkdir -p ${outputPath}/Kras_Final
thread=10
for tag in A B C D K ;do   # 1 2 5 10 M  W
    echo start; echo $Name
    Name=Sample_JZ24205452-NGS1105M4-${tag}
    FQName=$(echo "${Name}"|cut -d "_" -f2)

    FASTP=${outputPath}/fastp/${Name}
    ## fsatp clean
    echo "Q30_after: $Q30_after"
    # QC
    Q30_after=$(jq '.summary.after_filtering.q30_rate' ${FASTP}/${FQName}.json)
    # raw
    echo "raw and clean"
    FastqSize1=$(du -h ${samplePath}/${Name}/${FQName}_combined_R1.fastq.gz|awk '{print $1}')
    FastqSize2=$(du -h ${samplePath}/${Name}/${FQName}_combined_R2.fastq.gz|awk '{print $1}')
    
    raw_read_counts=$(fqtools count ${samplePath}/${Name}/${FQName}_combined_R1.fastq.gz)
    raw_base_counts=$(fqtools basetab ${samplePath}/${Name}/${FQName}_combined_R2.fastq.gz| awk '{sum += $2} END {print sum}')
    # clean  fastp
    clean_read_counts=$(fqtools count ${FASTP}/${FQName}_clean_R1.fastq.gz)
    clean_base_counts=$(fqtools basetab ${FASTP}/${FQName}_clean_R2.fastq.gz| awk '{sum += $2} END {print sum}')
    ############################### bwa 进行比对，使用前构建索引
    ## alignment fastp
    echo chr_distribution
    BAM=${outputPath}/BAM/${Name}
    ###  数据分析 
    # 染色体分布
    samtools idxstats ${BAM}/${FQName}.clean.umi.sort.bam > ${BAM}/${FQName}.chr_distribution.txt
    # 统计  去重前
    echo region and flagstat
    sambamba depth region -F "mapping_quality > 0 and not duplicate and not failed_quality_control" -t $(nproc) -L $BED ${BAM}/${FQName}.clean.umi.sort.bam  -o ${BAM}/${FQName}.c0.depth.txt
    sambamba flagstat -t $(nproc) ${BAM}/${FQName}.clean.umi.sort.bam > ${BAM}/${FQName}_base_flagstat.txt
    ## 上靶标率
    echo mapping_ratio
    mapping_ratio=$(cat ${BAM}/${FQName}_base_flagstat.txt | grep "mapped (" | awk -F '[()]' '{print $2}'|cut -d '%' -f1)
    ## 中靶率
    echo on target
    total_reads=$(samtools view -c ${BAM}/${FQName}.clean.umi.sort.bam)
    echo "total_reads: ${total_reads}"
    region_reads=$(samtools view -c ${BAM}/${FQName}.clean.umi.sort.bam -L $BED)
    echo "region_reads: ${region_reads}"
    on_target=$(echo "scale=4;($region_reads / $total_reads) * 100" | bc -l)
    ###  trim
    echo on trim
    MT65_reads=$(samtools view -c ${BAM}/${FQName}.lengthMT65.cat.bam)
    trim_percentage=$(( (total_reads - MT65_reads) * 100 / total_reads ))
    averageDepth=$(awk '{sum+=$5} END {print sum/NR}' ${BAM}/${FQName}.c0.depth.txt )
    ### dedup umi|gatk  去重后

    echo -e "sampleID\tfastq_size\traw_reads\traw_bases\tclean_reads\tclean_bases\tqc30\ttrim_percent(%)\tmapping_rate(%)\ton-target_percent(%)\tmean_depth\tmean_dedup_depth" >  ${BAM}/${FQName}_QC_output.txt
    echo -e "${FQName}\t${FastqSize1}/${FastqSize2}\t${raw_read_counts}\t${raw_base_counts}\t${clean_read_counts}\t${clean_base_counts}\t${Q30_after}\t${trim_percentage}\t${mapping_ratio}\t${on_target}\t${averageDepth}\t${deDupliAverageDepth}" >>  ${BAM}/${FQName}_QC_output.txt
    echo succeed

    mv ${BAM}/${FQName}_QC_output.txt ${outputPath}/Kras_Final
done
